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A method is proposed for the estimation of the transfer 
function of a linear, time-invariant system with no numera- 
tor dynamics, from random input and output data. The method 
employed utilizes the Fast Fourier Transform Algorithm and 
Least Squares Estimation to obtain the coefficients of the 
system's transfer function. 

The procedure has been modeled in FORTRAN IV on an 
IBM-360 computer. The results of simulation show the 
feasibility of estimating the order of the transfer function 
and its coefficients. 
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I . INTRODUCTION 



The problem of identifying a system (or plant) has been 
the object of much thought and effort as evidenced by the 
Identification in Automatic Control Systems Symposium [1] . 
The effort has been primarily expended on methods which are 
suitable to a time domain analysis. Identification in the 
frequency domain has in general involved lengthy computa- 
tions to obtain a Fourier Transform if the computations are 
done digitally, or a tradeoff between many parallel narrow 
band filters or repetitive playback of the system input 
and output if the analysis is done by analog methods. The 
method to be used in this investigation has been constrained 
to make use of the Fast Fourier Transform algorithm. There- 
fore digital computations and techniques logically follow. 
The Fast Fourier Transform (FFT) algorithm reduces the 
number of complex multiplications required in the evaluation 
of the Fourier Transform and thus may provide a means of 
performing system identification on-line or approximately in 
real time. The FFT is discussed in Appendix A. 



9 



II. PROBLEM STATEMENT 



The mathematical model of a linear, time-invariant 
system may be expressed as a differential equation or as 
a transfer function. The determination of the model 
parameters from observations on the system which is to be 
modeled is generally called an identification problem. 
This investigation considers the modeling of a system, 
which is characterized by a vector differential equation, 
such that the coefficients of the equivalent transfer 
function can be estimated. 

Let a linear system transfer function be written as 



X(f) = 1 . 

vtfT STfT 5 



The equivalent N 



th 



N ,n , , s 

d x j t _ ) . 

n=0 n dt n 



N 

D(f) = E d ( j a) ) n . (1) 

n=0 n 



order differential equation is written 



v(t ) . 



( 2 ) 



d N = 1.0 

From (2) the equivalent N-vector differential equation 
(3) can be obtained. 



; i 1 (t) ^0 10 

x 2 (t) 0 0 1 



Vl (t) 

x N (fc) 



0 

0 



0 0 0 



-d Q -d x -d 2 . . . -d N _ 1 



OR 



x 1 (t) 

x 2 (t ) 



X N-l (t) 

x N (t) 



v(t) 



( 3 ) 



x( t ) = Ax( t ) + Bv( t ) 
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The solution of (3) is 



t 

x(t) = $(t ,t 0 )x(t 0 ) + f $(t s x) • B • v( t ) dx j (4) 

t 0 

where x(t) is the system state vector and the element x-^(t) 
is the system output x(t). The discrete version of (4), 
which is derived in Appendix B, is used as the equivalent 

representation of (1) and the coefficients d^, d^, d^ 

are estimated. The estimation scheme is represented by 
Figure 1. 

The investigation is to be constrained so that v(t) is 
broad-band white noise and the Fast Fourier Transform 
algorithm must be used in the identification scheme. 
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UNKNOWN SYSTEM (EQUATION 4) 




Figure 1. Estimation Scheme Representation 
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III. SYSTEM MODEL FOR A DIGITAL COMPUTER 



A linear, time-invariant system represented by the 
vector differential equation (5) has equation (6) as its 
solution . 



x(t) = Ax(t) + Bv( t ) 



(5) 



x(t ) 



A(t-t n ) 
e x(t Q ) 




£ A(t 



( 6 ) 



If v(t) is deterministic and the convolution integral 

can be solved, then a discrete solution of (5) can be 

A( t-t Q ) 

modeled on a digital computer. e can be converted 

to a square matrix by using the Laplace transform method 
[2, pg. 316]. However, since v(x) is not deterministic the 
convolution integral cannot be solved, so an approximation 
to (6) is required. 

Considering only the unforced response of (6), 

A ( t-t n ) 

x(t) = e x(t Q ) . (7) 

At some instant of time t = t^ = tg + MT, 

— -( t M ) = £AMT — ^ t o ^ (8) 

which is identical to the unforced response of equation (3) 
of Appendix B where x(MT) = x(tg + MT) = x(tjyj). 
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( 9 ) 



AT 

x((m+l)T) = e x( mT ) 
x(T) = e AT x(0) = £ AT x(t Q ) 

x(2T) = e AT x(T) = e A2T x(tg) 

• • 

• ♦ 

• • 

x(MT) = e A ^^x(tg) (10) 

Therefore, since (10) is equivalent to (8) which is the 
discrete solution of (7), then (9) is the discrete solution 
for the unforced response of (6) for t = t Q + (m+l)T 

or in the notation of Appendix B, t = (m+l)T. 

The convolution integral from (6), 

g(t) = / e A ^ t-T ' ) Bv(T)dT, (11) 

t 0 

cannot be forced to yield an exact discrete solution without 
performing the indicated integration. Therefore, the sam- 
pled data approach in Appendix B is necessary. Normally 
this approach yields 

g((m+l)T) = / T e A(T " x) Bv(mT)dT, (12) 

0 

where v(mT) is assumed to be constant for mT<.t<(m+l)T 
[2, pg. 340]. This approach presupposes that the bandwidth 
of v(t) is small compared to Since v(t) is assumed to 

be broadband noise, T should be small. However, as T 
decreases, the number of iterations required in the simula- 
tion for a given time span of x(t) increases. Each iteration 

14 



requires that both the unforced and forced responses of 
equation (6) be computed. T should then be large to min- 
imize computation and small to allow the maximum transfer 
of information to the system concerning the variations of 
v(t). Equation (12) of Appendix B was chosen as a compro- 
mise. By using this equation, v(t) can be sampled at a 
rate greater than and the unforced response of equation 
(6) need only be computed once every T time interval. 
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IV. SYSTEM IDENTIFICATION SCHEME 



A continuous, time-invariant system is to be identified 
by estimating the coefficients of its transfer function. 

The input to the system is broad-band white noise of known 
mean and variance. The identification scheme is depicted 
by Figure 2, where F denotes the Fourier Transform, i.e., 

V(f) = F [v(t ) ] 

X(f) = F [x(t ) ] 

v(t) - broadband white noise 
x(t) - system output 

The object of this investigation is to make use of 
the Fast Fourier Transform (FFT, Appendix A) in place of F, 
and Least Squares Estimation in the computations for the 
transfer function coefficients. 




Figure 2. Identification Block Diagram 
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In Figure 2, v(t) is the output of a random process. 

In particular it is the output of a gaussian random process. 
Since v(t) is random, then x(t), V(f) and X(f) are random 
processes. V(f) and X(f) are not themselves deterministic; 
however, from Goldman, [ 3 , pg. 2791 

H ( f ) = X(f)/V(f) , (13) 

where H(f) is the frequency response of the unknown system. 
Since H(f) is assumed to be of the form 

H(f) = -N 1 DTfT 

E d ( j 2 7T f ) n 
n=0 n 

d N = 1 • 0 5 

then (13) can be written as [Appendix C, equation (2)] 

X(f ) • D(f ) = V(f) . (1*1) 

From (14) a residue r(f) is formed, 

r(f) = X ( f ) • D ( f ) - V(f). (15) 

Through the Fourier Transform, measurements on x(t) and 
v(t) yield X(f) and V(f). A modified least squares esti- 
mation procedure then uses r(f) from (15) and its complex 
conjugate to form 
F 

Q = / r(f) • r*(f)df (16) 

0 
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where r*(f) is the complex conjugate of r(f). The object 
is then to minimize Q with respect to the d 's, the 
coefficients of D(f). 

Since the PFT [Appendix A] is used in place of the 
Fourier Transform of Figure 2 and the unknown system is 
modeled in accordance with Appendix B, the identification 
scheme takes on the form of Figure 3. 

Since u(t) is broadband noise, the low-pass filter is 
a necessary addition to the identification scheme to prevent 
foldover of the alias frequencies of . The low-pass 
filter used is a discrete version of an R-C filter with 
its break point on a Bode plot being greater than the 
expected bandwidth of the system under investigation. The 
sampling rate of the inputs to the FFT ' s was then chosen 
to be at least twice as great as the low-pass filter cut- 
off frequency, that frequency at which the filter output 
drops by 3db . 

and VX are the results of the Discrete Fourier Trans- 
forms of x(t) and u(t) which are observed through a time 
window of t seconds, i.e. 

X. = X(f . ) = X(-) 

1 1 T 

V. = V(f . ) = V(-) , 

where t is the time interval over which u(t) and x(t) are 
observed. As the index i goes from -tW to xW, X^ and 
represent the complex values of the Fourier Transforms of 
x(t) and u(t) at the discrete frequencies -W to +W. 
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Figure 3, Identification Block Diagram with Signal Sampling 



The ' s and Vh ' s exist only at frequencies which are 

multiples of If x(t) and v(t) contain only frequencies 

which are exact multiples of ^ then the ' s and ' s 

will be the exact Fourier coefficients of x(t) and v(t)[6]. 

However, since it is unlikely that the only frequencies 

present are multiples of these other frequency components 

will cause the X.'s and V.’s to be altered. It would be 

1 1 

advisable to weight the frequency spectra of x(t) and 
v(t) with the hamming or hanning functions [4, 5, 6, Appen- 
dix A] to minimize the effects of these other frequency 
components. This weighting was not performed, although it 
is admittedly advisable, since computer computation time 
would be increased. The trade-off of less computation time 
for increased accuracy was felt to be justified since much 
of the computer time used in this investigation had to be 
charged to other projects. 
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V. RESULTS 



Many of the variables in the computer simulation of the 
"unknown" plant and the identification scheme [Table Dl, 
Appendix D] are peculiar to the approach taken in this 
investigation. For an actual application of the identifi- 
cation scheme the "unknown" plant is not simulated and only 
the sampling rate and number of samples to be taken must 
be determined. The expected or known signal bandwidth 
dictates the sampling rate, where the minimum rate is the 
Nyquist sampling rate. An increase in the number of 
samples taken increases the frequency resolution of the 
FFT output and the accuracy of the estimated transfer 
function coefficients. However, it also increases compu- 
tation time and computer memory required. 

A tenth order system such as equation (17) was assumed 
for each trial. 



X(S) _ 1.0 

vTsT 10 „ 

E d S 



(17) 



The low-pass filter shown in Figure 3 is represented by 
equation ( 18 ) . 

V(S) _ 2tt ' 100 nfn 

uTsT S+2 tt • 100 K 1 

Table 1 presents the most encouraging of the results 
of this investigation. Trials 1 and 2 represent first-order 
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systems with real S-plane poles. Trials 3 and 4 represent 
second-order systems with real and complex poles respect- 
ively. Trial 5 represents a second-order system with 
complex poles and trial 6 is a third-order system with 
real poles of - 10 , -20 and - 30 . Only the five lowest-order 
coefficients of the estimates are given in Table 1 since 
the higher-order estimates become negligibly small. The 
results given are the averages of fifty consecutive esti- 
mates made for each system's coefficients. 
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TABLE 1 



ESTIMATION RESULTS 
Actual System coefficients = d n 

A 

Estimated coefficients = d^ 



d^ through dg = 0 



TRIAL 


#SAMPLE 

POINTS 


SAMPLING 
RATE (HERTZ) 


n 


d 

n 


A 

d 

n 

LOW PASS FILTER 
IN OUT 


1 


256 


166 


7 


0.0 


9.6 -10" 9 


-2.6*10-8 








8 


0.0 


4.2-10" 6 


1.0-10" 5 








9 


0.0 


-1.8*10"3 


-3.5-10 -3 












_1 


_1 








10 


1.0 


9.1*10 


9.3-10 


i 




■ 


11 


6.0-10 1 


6.2-10 1 


6.2-10 1 


2 


i 

512 


416 


7 


0.0 


3 

O 

i — 1 
i — 1 

i — 1 

1 


-9.2-10 -9 








8 


0.0 


2.0-10 7 


1.2-10" 6 








9 


0.0 


-4 

-4.2-10 


oo 

o 
1 — 1 

VO 

C\j 

1 








10 


1.0 


7.9* 10 -1 


1 — 1 

1 

o 
1 — 1 

VO 

l n 








11 


6.0-10 2 


6.1-10 2 


6.8* 10 2 


3 


512 


250 


7 


0.0 


4 .7-10 -6 


5.3-10 -6 








8 


0.0 


-4 

-9*3-10 


-8.7-10 -4 








9 


1.0 


i — 1 

1 

O 

i — 1 

VO 

CO 


8.6 -10 _1 








10 


8 . 0 ‘ 10 1 


7.5-10 1 


7.6 -10 1 








11 


1.5-10 3 


1.5-10 3 


1 .5 -10 3 


4 


512 


83 


7 


0.0 


5.3-10" 6 


2.3-10 -5 








8 


0.0 


-4 

-1.3-10 


-1.5-10 -3 








9 


1.0 


9. 3-10" 1 


1.0 








10 


6.0-10 1 


5.0-10 1 


5 .1 -10 1 








11 


on 

O 

i — 1 

VO 

on 


3 . 6 • 10 3 


3 . 7 " 10 3 
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TRIAL 


^SAMPLE 

POINTS 


SAMPLING 

RATE(HEFTZ) 


n d 

n 


A 

d 

n 

LOW PASS FILTER 
IN OUT 


5 


512 


4l6 


7 0.0 

8 0.0 

9 1.0 

10 1.2-10 1 

11 3.6-10 3 


2.1-10 -6 
-5. 3‘10 -5 
6 . 2 ‘ 10 -1 
1.1-10 1 
2 . 4 ' 10 3 


1.3*10 6 

-4 

-1.5*10 

5.4-10" 1 

8.2 

2 . 1 • 10 3 


6 


1024 


83 


7 0.0 
8 1.0 
9 6.0-10 1 

10 l.l'lO 3 

11 6.0-10 3 


2.5* 10 -3 
7.8* 10 -1 
5 . 8 * 10 1 
1.0-10 3 
6.3*10 3 


2.0-10 -3 
7.8* 10 -1 
5.6-10 1 
1.0-10 3 
6.1-10 3 
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VI. CONCLUSIONS 



The identification scheme presented in this investiga- 
tion does provide good estimates of the transfer function 
coefficients for an unknown system with no numerator 
dynamics. The estimated coefficients in Table 1 are seen 
to be within 50% of the actual coefficients for all cases 
except when the actual coefficient equals zero. For these 
coefficients the estimates rapidly approach zero since the 
estimates are for higher order coefficients than those that 
exist in the actual system. The identification scheme 
therefore provides an indication of the order of the system. 

The need for a low pass filter on the input to the 
"unknown" system has not been substantiated. Trials 1 
through 6 of Table 1 indicate that the filter can be deleted 
with no disastrous effects. 

The feasibility of using the Fast Fourier Transform as 
a tool for system (plant) identification has been shown. 
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VII. RECOMMENDATIONS 



This investigation has been unduly long and involved 
due to the requirement that all the modeling be done on a 
digital computer. Many of the parameters listed in Table 
D1 of Appendix D required a large number of trial runs to 
determine a best set of values for the overall simulation. 
These trials were necessary to aid in finding the minimum 
iteration rate and input signal sampling rate required to 
adequately simulate the continuous "unknown" plant since 
these rates are fundamental to the total simulation time 
(not to be confused with identification time). It is 
suggested that any further investigation in this area be 
done on a Hybrid computer so that the continuous plant may 
be run on the analog portion of the computer and the iden- 
tification scheme on the digital portion. 

As a starting point the identification scheme requires 
that a guess be made as to the maximum number of coeffi- 
cients in the unknown plant's transfer function. A major 
effect of this guess is that it determines the dimensions 
of matrix W of equation (15), Appendix C. Since the inverse 
of W is required it would appear that its dimensions should 
be kept as small as possible. However, it has been observed 
that the dimensions of W that produce the most nearly cor- 
rect estimates of the transfer function coefficients are 
not the smallest (consistent with the number of coefficients 
known to exist) nor the largest (which depend on the amount 
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of computer memory reserved for W) . The best dimensions 
for W are different even for different plants of the same 



order. It is suggested that any future 
izing this identification scheme should 
of avoiding the inverse of W or attempt 
of W -1 . 



investigation util- 
involve a method 
an analytic solution 
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APPENDIX A 



The Fast Fourier Transform (FFT) 



The FFT is used to evaluate the complex Fourier trans- 
form of a sequence of complex numbers. It is a finitie 
discrete Fourier transform (DFT) with the restriction that 
the number of samples in the complex sequence is constrained 
to be an exact power of two. The power-of-two constraint 
on the number of samples allows implementation of an algor- 
ithm on a digital computer for rapid computation of the 
finite DFT. 

The mechanics of the FFT and its characteristics has 
been widely discussed in the literature [7> 8, 9 , 10, 11, 

12, 13]. A few of the characteristics will be presented 
here . 

From Goldman [3] , a function g(t) which exists only for 
— T T 

^- < t <— and which has a Fourier transform G(f) which exists 
only for -W <f <W can be expressed by the DFT pair 




( 1 ) 




. 2 mn 
J 2TW 



( 2 ) 



where g(^j) = g(t), for t = ^ 



g(t) and G(f) are completely determined by their 2TW+1 

sample values g(^) and G(^) respectively. 

The assumptions leading to (1) and (2), that g(t) exist 
-T T 

only for — < t < ^ and G(f ) exists only for -W < f < W, are 
physically unrealizable. For example, consider the case 
where g(t) is a section of a continuous function x(t) such 
that 

x(t) is continuous for < t < «> 
g(t) = x ( t ) • y ( t ) 

y ( t ) = u(t+|) - u(t-|) 
u(t) is the unit step function . 

Then G(f) = f°° X(q)Y(f-q)dq 

— oo 

and Y( f ) = T sln TI ^ fT ) . 

Y(f) and therefore G(f) are continuous for -°° < f < 00 . 

By observing (sampling) a section of x(t) such that 
-T T 

g(t) = x(t) for — <t <— the assumption leading to (1) 
and (2) are violated. 

If a y(t) in (3) is chosen so that g(t) = x(t)* y(t) 

-T T 

has negligible value outside of the range — < t < ^ and G(f) 
is negligible outside the range -W < f < W, then (1) and (2) 
can be made to be approximately exact. y(t)'s to satisfy 
the above include the Hanning and Hamming functions [4, 51 
and the Parzen spectral window [14]. 



(3) 

(4) 

(5) 
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For the FFT, (1) and (2) are normally expressed as 



g(m) 

G(n) 



N-l j 

S G (n ) e 
n=0 



27Tmn 

N 



N-l 

N £ n « <m) 

m=0 




2irmn 

N 



( 6 ) 

(7) 



where g(m) = g(^) 

G(n) = G(^-) 

n,m = 0 to N-l 
N = 2 k . 

Some authors put the weighting j-j- on (6) rather than (7). 
However, the FFT used in this investigation is from the 
IBM/360 Scientific Subroutine package ( 36 OA-CM-O 3X) and in 
this algorithm the jjj- weighting is placed on (7). 
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APPENDIX B 



Discrete Solution for x(t) 

Given the vector differential equation (1) it is desired 
to obtain a discrete solution for x_(t). 

x(t) = Ax(t) + Bv(t) j (1) 

T 

where x(t) = [x^t) x 2 (t) x N (t)] 

T 

x(t) = [x 1 (t) x 2 (t) x N (t)] 

A is an (N,N) constant matrix 
B is an (N,l) constant vector 
u(t) is a scalar . 



The superscript T denotes the transpose of the vector. From 
Ogata [2] the continuous time solution of (1) is 



x(t ) 



A(t-t Q ) 

e 



t 

x(t n ) + / 

u t 
r 0 



e A(t-T) 



Bv( T ) 



dT, 



(2) 



A(t-t Q ) 

where e is the transition matrix $(t,tQ) 

For a discrete solution of (2), where x(t) is sampled 
at intervals of time T, let 



t 



0 



nT 



t = (n + 1)T . 



Then 



x [(n+l)T] = e AT x(nT) + / ^ £ A[(n+l)T t] 

nT 



Bv( x ) dx (3) 
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The usual procedure Is to assume that v(x) = v(nT) for, 
nT <. x <(n + l)T. The implementation of v(x) in this inves- 
tigation possesses characteristics which make it advisable 
to sample this forcing function at as high a rate as 
possible. Since v(x) is assumed to be broadband noise, it 
is varying very rapidly with respect to time. Therefore T 
should be very small if v(x) is to be fairly constant over 
the sampling interval. Also it must be noted that in the 
simulation of the identification problem v(x) is represented 
by a sequence of numbers from a random number generator. 

The magnitude of the numbers in this sequence are indepen- 
dent of T in that if five numbers are produced they will 
have the same magnitudes whether they are produced at 
intervals of T or T/2 . The implication then is that for 
any two consecutive numbers from the generator, such as v(i) 
and v(i+l), those which are produced at T/2 second intervals 
represent a v(x) that has a higher rate of change with 
respect to time (and a higher frequency) than those that 
are produced at T second intervals, i.e., 

v(i+l)-v(i) ^ v(i+l)-v(i) ' 

, T/2 T 

Therefore if the v's in (3) can be generated at intervals 
smaller than T, then the effective bandwidth of the noise 
can be made wider as the sampling interval of v(x) becomes 
smaller. Some averaging of the rapid fluctuations of v(x), 
weighted by the system parameters in the matrix A, will be 
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made by the integral in (3). The result is that v(t) can 
be sampled at a high rate (higher than ^-) and the unforced 
response in (3) need only be computed every T seconds. To 
implement this approach let. 



x [ (n+1 )T] = 



(m+1) 

flT M-l M 

e R x(nT) + Z / 

m=0 (^n)T 



+n)T 

[ (n+1 )T-t ] B v( t ) dx 



(4) 



v(t) = v[(|J|+n)T], (^+n)T < t <( m ^ 1 +n)T . 

Since B is a constant vector and v(x) will be a constant 
over the interval of integration. 



x [ (n+1 )T] 



AT / m x , 
e x(nT) + 



M-l 

Z 

m=0 



(2^H-n)T 

/ eA[( n+ i) T _x ]dx-B-v [(jJ+n)T] 



(^ n )T 



(5) 



To remove the parameter n from the integral in (5) let 



x = nT-A 
A = nT-x 



Then 



_( d 1+ ^ )T 

M-l ^ M ; 

x[(n+l)T]=e AT x(nT)+e AT Z -/ e AX dA • B • v [ ({Jj+n ) T] 

m=0 -mT 
M 



(6) 



To remove the negative limits on the integral in (6) let 



A 

x 



mT 
T M 

, mT 
A M 
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Then (6) becomes 



am AT M-l M -A(t+^) 

x[(n+l)T] = £ Ai x(nT)+e fl E / e dx *B • v [ (^fn)T] . 

m=0 0 



(7) 



Since the integral in (7) is not a function of m, then 



M M-l m 

x[(n+l)T] = e AT x(nT)+e AT / e _AT dT E e~ A ^Bv[(^+N)T]. (8) 

0 m=0 



To evaluate the integral in (8) on a digital computer. 



—At 

e is expanded in an infinite series 



T 

M 



M 00 



r -At, v (-At) 1 , 

/ e dT = / E - — - dT 



0 



0 i=0 



i! 



T 

00 ( -A ) i , M i 

= E - — — * / T dT 

i=0 1 ‘ 0 



= E 

1 = 0 



( -A ) ' 
i! 



T (i+l) 
M 



(i+l) 



.-AT/ 
T ” M 
M l (i+l) ! 



T 

M 



-At. 



Equation (9) is the solution of / e dT to be used, 

0 

ever, if A is known to be nonsingular, then 



(9) 



How- 



T 

M 

f -At , „-l rT _-AT, 

J e dT = A [I-e-p— ] 

0 



( 10 ) 



where I is the identity matrix of the same dimensions as A, 
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A - "*" should exist except for those cases where the system 
characterized by (1) has a pole at S=0. However, to avoid 
compounding errors by taking an unnecessary matrix inverse 
on the computer, (9) is the solution of the integral to be 
used . 

Prom (8) and (9) and letting A = j^. 



x [ (n+1 )T] = 



AT / rn \ i AT ~ 

e x (nT ) +e • A • E 

i = 0 



OR 



( -AA ) 1 M t , 1 -AmA „ , m , , N 

, • £ e *B*v(nT+mA) 

u ; • m=0 

( 11 ) 



where 



x [ (n+1 )T] = $(T) x(nT) + T (T ,M) v(nT ,M) 



( 12 ) 



x(nT) is the (N,l) state vector of the plant at 
time t=nT 

AT 

$(T) is the (N,N) transition matrix e 
T is the transition time between the system's states 



v(nT,M) is the (M,l) control vector 



v(nT) 
v(nT+A ) 



v(nT+(M-l) A) 



f(.T,M) is an (N,M) transfer matrix such that 



f (T,M) = $(T) 



A 2 



(-AA) 



i = 0 ( i+1 ^ 



-AA 



-2AA 



.-(M-l)AA . B 
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APPENDIX C 



Implementation of the System Identification Solution 



A given system is represented by Figure Cl, where 
K k 

DO) = I d, (jco) • It is desired that the coefficients 
k = 0 

of DO) be identified through measurements on v(t) and x(t). 



v(t ) 




1 


x(t ) 


VO) 




DO) 


xO) 



Figure Cl. System Representation in the 
Frequency Domain 



Since 



1 _ XO) 

dTToT vOT ’ 



(i) 



then 



XO) DO) = VO) 



(2) 



If XO) and VO) are known for discrete values of eo where 
w -*■ to. then 

l 

X(o) i )D(o) i ) = V0 1 ) • (3) 



XCoi. ) andVCco^) result from application of the Fast 
Fourier Transform algorithm [Appendix A] to measurements 
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of x(t) and v(t). Therefore X(oo^) and V(un) are assumed 

A 

to be known. If D(uk) in (3) is replaced by D(uk), the 
estimate of D(uk), then (3) becomes 



X(w.)D(to.) = V(w i ) (4) 



Let r(to.) equal the residue, or error, due to the 

/v 

estimate D(to^ ) . 



r(to.) = X(to.)D(w.) - V(u>.) (5) 



Then applying a least squares estimation procedure to 



I 

Q = Z r (to. ) • r* (to . ) , (6) 

1 = 0 1 1 

* denotes the complex conjugate 
Q is to be minimized. 

From (5 ) and (6 ) , 

Q= E [X(w. )D(w. )-V(oj. ) ] • [X ( to . ) D ( to . ) - V ( to . ) ] * 

i = 0 x x x xxx 

(7) 

Let 



X(oo. ) = X. 

l l 



V(co 1 ) = V. 



D(w j _) = Z d k (jto 1 ) k = D 1 
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Then, setting the partial of Q with respect to a coeffi 

cient d equal to zero, 
m ’ 



90 I . 3D* 3D. 

= Z (X D -V . )X* + X (X D -V. )* = 0, 

3d i=0 3d 3d 11 

m mm 



where 



3D. 

"T 1 = (J^) 



m 



3d 



m 



3D*. . K A , 

1 3 - 5 / • \k / , N.rri/ . N m 

— — = -a- 2 d.(-ju. ) = (-1) (ju.) 

3d 3d k=0 x 

m m 



From (8), ( 9 ) and (10), 



I ^ 

l X ± X*. [D 1 (-l) m (jw i ) m + D* (jcd,) 1 "] 
i = 0 

I 

Z V.X*. (-l) m (joj. ) m + X.V*. (ju>. ) m . 

. A l l ° l l l ° l 

i=0 

Since (j ) m is common to both sides of (11) and is not a 
function of i. 



Z X.X*. [ Z (-l) m d, (jw. ) k + (-l) k d, ( jw. ) k ] w . 

. _ i i . A k d i ' k 0 i i 

1=0 k=0 

I 

= Z [(-l) m - X*.V. + V*.X.] co. m , 

ii il l 

1 = 0 



( 8 ) 



( 9 ) 



( 10 ) 



(ID 



( 12 ) 
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where has been replaced by its summation. The ' s are 
not functions of i or m, and 

(-l) m +(-l)^ = 0, if m and k are not both even or odd 

= 2, if m and k are both even 

= -2 if m and k are both odd 



Therefore, for m = 0 to M and k = 0 to K the left side of 
(12) can be expressed as 




0 

w(l,l) 

0 

w(3,l) 



w (M , 1 ) 



w( 0 ,K ) 
w(l,K) 



• d (13) 



w(M,K)_ 



where 



w 



(m,k) = j k [ (-1 ) m +(-l ) k ] 



d = 



d 



K 
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The right hand side of (12) can be expressed as 



i = 0 



Re (X.*V.) 

1 1 

-jIm(X.»V. ) . 

o I. 1 l 

Re (X.*V.) 

1 1 1 



(- 1 ) M -X. *v.+v.+v. *x. • . M 

11111 1 



( 1*0 



where Re(X.*V.) = real part of X.*V. 

ii ^ ii 

Im(X^*V^) = imaginary part of X^*-V^ 

co . = 2 Trf . = ^ 11 . 
i it 

t = the time span over which v(t) and x(t) are 
observed . 

Now to go from the foregoing general derivation to the 
particular case where the order of the system is N, let 
M = K = N, Then from (12), (13) and (1*1) and after simpli- 

A 

fication, the solution for d becomes 



~ -1 -] 
d = G W Z 



(15) 



where 



d = 



N 
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g 



0 



g l 




g 2 




g 0 = 1 



g i = J? 



g n g n-2 ^ 2 tt ^ 



1 
0 

w ( 2 j 0 ) 

0 

w( 4 ,0 ) 
w(N j 0 ) 

1=0 

w(m,k) = 0 j if m+k Is odd 

.m+k ,, 
l , otherwise 




0 w(0 ,N) 

w(l,l) w(l,N) 

0 - 

w(3,l) - 

0 _ 



w (N j 1 ) w(N,N)_ 



m,k = 0 to N 
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I 




i=0 



Re(X^*V^)i n for n even 

z = 
n 

-Im(X.*V.)i n for n odd 
1 1 

n = 0 to N 



The implementation of the identification scheme then 
takes on the form of Figure C2. 

For example, if the order of the system N equals 2 and 
v(t) and x(t) are observed through a time window of x 
seconds then. 




d 



1 



1 

0 



2 



0 



0 



0 



0 -(==■) 




1 

0 

. 2 

l 



0. 

. 2 

l 

0 




Re (X . *V. ) 

1 1 

-Im(X*V. ) i 
1 1 

Re (X . *V. )i 2 

_ l l _ 
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Figure C2 . Implementation of 
the Identification Scheme 



ANALOG 



DIGITAL 



d 



^3 



APPENDIX D 



Computer Program PLEST (PLANT ESTIMATION) 

The computer program -is composed of the main program 
and subroutines GAUSS, DHARM, PLANT, RECIP and EA . GAUSS 
and DHARM are called from the IBM/ 36 O Scientific Subroutine 
package ( 360A-CM-03X) and are not included in the program 
listing . 

Subroutine GAUSS is used to obtain a sequence of real 
random numbers with a normal distribution and a selectable 
mean and standard deviation. The expected value of the 
power spectrum of the output of GAUSS has been investigated 
by utilizing the Fast Fourier transform and has been found 
to approach a constant value over all frequencies. GAUSS 
is therefore used to approximate a source of "white" noise. 

Subroutine DHARM computes the Discrete Complex Fourier 
Transform of a sequence of numbers. DHARM is used to 
obtain the frequency spectrum of the input and output of 
the system (plant) for which the transfer function coef- 
ficients are to be determined. 

Subroutine RECIP is used to obtain the inverse of a 
square matrix with dimensions of from (1,1) to(12,12). 

Subroutine EA computes either. 
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( 1 ) 



or 




( 2 ) 



where A is a square matrix 

When the sum of the magnitudes of all of the elements 



Subroutine PLANT computes the 0 and r matrices of 
equation (12) in Appendix B. The inputs to PLANT consist 
of either the coefficients of a transfer function for a 
system or the A matrix and B vector from a vector differ- 
ential equation such as; 



The MAIN program (PLEST) performs the normal bookkeeping 
functions required of any comprehensive program. It also 
utilizes subroutine PLANT so that a continuous system such 
as Figure D1 is converted to the discrete system in 
Figure D2 . The MAIN program, subroutines DHARM and RECIP 
perform the operations required by Appendix C to identify 
the system under test. 

The primary program parameters are identified in the 
table of Computer Program Data Cards (Table Dl). 



in the N^ h term of (1) or (2) is less than or equal to 
-20 

10 , the computation is stopped. 



x(t) = Ax(t)+Bv(t) 



( 3 ) 
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TO FOURIER 
TRANSFORM 



TO FOURIER 
TRANSFORM 



Figure D1 . Continuous Version of 
the System Under Test 



H6 



u(T- n) 



NOISE GENERATOR 
(FROM SUBROUTINE 
GAUSS) 






DELAY 




r 




(T-N) 












$ ( T • N ) 


^ . 









N u(T-(N-l)) 



LOW PASS FILTER 



one cy cTe) 

v(T*N *n) 



l 




x(T-N-M-n) 



M v(T-N-(N-0)) SYSTEM UNDER TEST 



l 



T • N • M • K 



^/T-N-M- 



TO FAST 
FOURIER 
TRANSFORM 



L-l 

NOTE: I r ( !•£) - y(r'l) 

11 = 0 

= P( x (nx ,L) 

from Equation (12) 

Appendix B 

TO FAST 
FOURIER 
TRANSFORM 



K 



l 



Figure D2 . Discrete Version of the System under Test 
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TABLE D1 



COMPUTER PROGRAM DATA CARDS 



Card Format Name 



Comments 



1 14 NJOBS 

2 10A8 

3 DIO. 4 YARN 



4 DIO. 4 T 

5 14 N 



6 1 4 M 



7 1 4 K 



The number of sets of data cards 
to be operated on. 

Arbitrary comments to be read and 
written on the print-out. 

The variance of the random signal 
generated by the subroutine gauss 
(the mean has been preset to 0). 

The time (T) between samples for 
the filter input. 

The number of successive inputs 
to the filter for one output. 

The iteration time of the filter 
(TN) is determined by TN=T.N. 

0 < N < 12 

The number of successive inputs 
to the plant under test for one 
output. The iteration time of 
the plant is determined by 
TMN = T-M-N. 0 < M < 12 

For every iteration of the 

plant under test the plant input 
and output are sampled once for 
the inputs to the FFT . The inputs 
to the FFT occur at intervals of 
TKMN = T • K • M • N . 0<K<12 
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Card Format Name 


Comments 


8 I 4 NRUNS 


The number of times the plant 
coefficients will be deter- 
mined. The resulting coeffi- 
cients are the average of 
the results from NRUNS esti- 
mates of the coefficients. 


9 14 LTRMS 


A lower-limit guess on the 
number of terms in the plant 
transfer function. 

1 < LTRMS < NTRMS 


10 14 NTRMS 


An upper-limit guess on the 
number of terms in the plant 
transfer function. 

1 < NTRMS < 12 


11 I 4 ISTST 


The number of iterations the 
filter-plant combination will 
go through to allow it to 
reach a steady state condition 


12 14 M ( 1 ) 


2^^ samples of the plant 
input and output are to be 
used on each run (each NRUNS 
from card 8) to estimate the 
transfer function coefficients. 
3 £ M ( 1 ) < 10 


13 12 IFLT 


The method by which the fil- 
ter will be characterized.. If 
IFLT = 1 the transfer function 
coefficients are to be entered. 
If IFLT =2 the bottom row of 
the A matrix and each element 
of the B vector is to be 
entered . 
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Card Format 



Name 



Comments 



1^ 12 NDPFLT The exponent of highest power 

of S in the denominator of 
the filter transfer function. 
0 £ NDPFLT £ 11 

15 12 NPFLT The exponent of highest power 

of S in the numerator of the 
filter transfer function. 

0 < NPFLT < NDPFLT 



NOTE: if IFLT = 2, go to (2) 16A 



(1) 


16A 


DIO. 4 


DE ( 1 ) 


(1) 


16b 


DIO. 4 

• 


DE ( 2 ) 

♦ 


(1) 


16- 


• 

DIO. 4 


• 

DE ( NDPFLT+1 ) 


Cl) 


17A 


DIO . 4 


UN ( 1 ) 


(1) 


17B 


DIO. 4 

• 


UN ( 2 ) 


(1) 


17- 


• 

DIO . 4 


UN (NPFLT+1 ) 


NOTE 




if IFLT 


= 1, go to 18 


(2) 


16a 


DIO. 4 


A (NDPFLT, 1) 



The coefficient of the highest 
power of S in the denominator. 



The coefficient of in the 
transfer function denominator. 

The coefficient of the highest 
power of S in the transfer 
function numerator. 



The coefficient of in the 
transfer function numerator. 



The element A(NDPFLT,1) in 
the filter A matrix. 
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Card 


Format 


Name 


Comments 


(2) 


16B 


D10. 4 


A( NDPFLT ,2) 




(2) 


16 - 


D10. 4 


A ( NDPFLT , NDPFLT ) 


The element A ( NDPFLT,NDPFLT) 
in the filter A matrix. 


(2) 


17A 


^r 

o 

i — l 

Q 


B ( 1 ) 


The element B(l) of the 
filter B vector 


(2) 


17B 


D10. 4 

• 


B ( 2 ) 




(2) 


• 

17- 


D10. 4 


B( NDPFLT) 


The element B(NDPFLT) of 
the filter B vector. 




18 


12 


IPLT 


The method by which the 



plant will be character- 
ized. Same comments as 
for card 13- 



Cards 19 through 22 (for the plant) are similar to 
cards 14 through 17 (for the filter). 



NOTE: to omit the filter enter the following data 



Card 

5 

13 

14 

15 

(1) 16A 
(1) 17A 



Data 

0001 

01 

00 

00 

1 . 0000D+00 
1 . 0000D+00 



Repeat Data Cards 2 through 22 NJOBS (from Card 1) times. 
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